Two Snuke
問題
https://gyazo.com/f9dd32211e6deb0099516ea09a18ca4f
問題条件でNが与えられる
N以下のnについて、それを5つの数に分ける
この一つに注目する、数をmとする
このmを2つの数i, jに分ける
$ i + j = m
ここで $ 0 \le i < jという制約がつけられている
とりうるすべてのi,jについての差の和を求めたい
$ f(m) = \sum_{i,j} j - i
ここまでは自力で求めたのだが、ここから先がわからなくなった
code:python
def f(x):
return (x + x % 2) * ((x + 2) // 2) // 2
形式的べき級数で解く
この関数fは自然数を引数に取る関数なので、数列と解釈できる
$ F = \sum_{m=0}^\infty f(m)
僕のfの実装では剰余や切り捨てが発生していたが、これを有理式で表すことができれば形式的べき級数の道具が使える
数列を観察
code:python
array([ 0, 1, 2, 4, 6, 9, 12, 16, 20, 25, 30,
36, 42, 49, 56, 64, 72, 81, 90, 100, 110, 121,
132, 144, 156, 169, 182, 196, 210, 225, 240, 256, 272,
289, 306, 324, 342, 361, 380, 400, 420, 441, 462, 484,
506, 529, 552, 576, 600, 625, 650, 676, 702, 729, 756,
784, 812, 841, 870, 900, 930, 961, 992, 1024, 1056, 1089,
1122, 1156, 1190, 1225, 1260, 1296, 1332, 1369, 1406, 1444, 1482,
1521, 1560, 1600, 1640, 1681, 1722, 1764, 1806, 1849, 1892, 1936,
1980, 2025, 2070, 2116, 2162, 2209, 2256, 2304, 2352, 2401, 2450,
2500])
僕は一般項をコードで書いてるので偶奇が影響するのわかってるから分けてみる
code:python
array([ 0, 2, 6, 12, 20, 30, 42, 56, 72, 90, 110,
132, 156, 182, 210, 240, 272, 306, 342, 380, 420, 462,
506, 552, 600, 650, 702, 756, 812, 870, 930, 992, 1056,
1122, 1190, 1260, 1332, 1406, 1482, 1560, 1640, 1722, 1806, 1892,
1980, 2070, 2162, 2256, 2352, 2450])
array([ 1, 4, 9, 16, 25, 36, 49, 64, 81, 100, 121,
144, 169, 196, 225, 256, 289, 324, 361, 400, 441, 484,
529, 576, 625, 676, 729, 784, 841, 900, 961, 1024, 1089,
1156, 1225, 1296, 1369, 1444, 1521, 1600, 1681, 1764, 1849, 1936,
2025, 2116, 2209, 2304, 2401, 2500])
この観察の時点でmaspyさんは$ F = \frac{P}{(1-x^2)^3}が確定したと書いているが僕にはそれがわからなかったので噛み砕く
一つ飛ばしの数列
1, 0, 1, 0, ... という数列を形式的べき級数でどう表現するのか
素朴に式を書いてみると
$ G = 1 + x^2 + x^4 + \ldots = \sum_{n=0}^\infty x^{2n}
$ G = 1 + x^2 + x^4 + \ldots
$ x^2G = x^2 + x^4 + \ldots
$ (1-x^2)G = 1
$ G = 1/(1-x^2)
一つ飛ばしの平方数
0, 1, 0, 4, 0, 9, ... という数列を形式的べき級数でどう表現するのか
これは$ H = \sum_{n=0}^\infty ((n+1) / 2)^2 x^nの奇数番目だけ取ったやつ
なお数列が1-originか0-originかで混乱したので0に揃えました
引いてみた $ F - H
code:python
array([-0.25, 0. , -0.25, 0. , -0.25, 0. , -0.25, 0. , -0.25,
0. , -0.25, 0. , -0.25, 0. , -0.25, 0. , -0.25, 0. ,...
つまり$ F = G / 4 + H
計算が面倒なので全部4倍して4H=Jとしておく
$ 4F = G + J
$ J = \sum_{n=0}^\infty (n+1)^2 x^n
初回は計算ミス をしててこの先で正しくない式を出してしまっていた code:python
array([ 1, 4, 9, 16, 25, 36, 49, 64, 81, 100, 121, 144, 169,
196, 225, 256, 289, 324, 361, 400])
array([ 1, 3, 5, 7, 9, 11, 13, 15, 17, 19, 21, 23, 25, 27, 29, 31, 33,
35, 37, 39])
Out37: array(1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2) Out38: array(1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0) $ J = \frac{1 + x}{(1-x)^3}
(1-x)^3は3回差を取ることに対応していたのか
$ 4F = G + J = \frac{1}{1-x^2} + \frac{1 + x}{(1-x)^3} = \frac{(1 - x)^2 + (1 + x)^2}{(1 + x)(1 - x)^3} = \frac{4x}{(1 + x)(1 - x)^3}
もっと簡単な式を求めて
maspyさんのやってる「分母に(1-x)が来そうだから掛けてみる」を試す 数列を有理式にする 形式的べき級数の積は、数列の畳み込み
np.convolveを使えば1行で計算できる
code:python
array([ 0, 1, 0, 0, 0, 0, 0, 0, 0, ...
(1 - x)(1 - x)(1 - x)(1 + x) を掛けたらxが残った
$ \frac{x}{(1-x)^3(1+x)}
僕が途中で計算間違いをしてたのだろうか？後で検証する。
ちなみにmaspyさんの当初の話どおり$ (1-x^2)^3を掛けてみたらこうなった
code:python
In 121: reduce(np.convolve, 1, 0, -1], 1, 0, -1, [1, 0, -1, fs) array([ 0, 1, 2, 1, 0, 0, 0, 0, 0, ...
$ \frac{x(1 + 2 x + x^2)}{(1-x^2)^3} = \frac{x(1 + x)^2}{(1-x)^3(1+x)^3} = \frac{x}{(1-x)^3(1 + x)}
整理すれば結果は同じ
Fが求まったとする
$ F = \frac{x}{(1-x)^3(1+x)}
Fは5つに分けたブロックの1つについてのべき級数だった
これはmを引数とする関数と解釈できて、合計がmである時の答えだった
$ [x^m]F = f(m) = \sum_{i,j} j - i
$ [x^n]F^5 が、5個のブロックの総和がnの時の答え
問題条件は「総和がN以下の時」なので、得たい最終的な答えは
$ [x^n] \sum_{n=0}^{\infty} F^5
この和を取るところと(1-x)で割ることはどう関係するのだっけ？？？
(1-x)で割るということは$ \sum_{n=0}^\infty x^n を掛けるということ
形式的べき級数の積は、数列の畳み込み。すべて1の数列と畳み込んでいるので、畳み込み結果の第N項は、元の数列のNまでの和になる
$ \sum_{n=0}^{\infty} F^5 = F^5 / (1-x) = \frac{x^5}{(1-x)^{16}(1+x)^5}
求めるべき形式的べき級数が決まったので、次はそれを解く
まず形式的べき級数の偶数次元/奇数次元だけ通す関数を定義する
$ \mathcal{E}(G(x)) = \frac{G(x) + G(-x)}{2}
$ \mathcal{O}(G(x)) = \frac{G(x) - G(-x)}{2}
要するにx軸で鏡映して足し合わせると奇関数は打ち消しあって消滅するということ
ここは説明のために(x)を明記したけど以下では省く, $ G^* = G(-x) とする
$ [x^n]F = P/Q を求めたい(Q(0)=1)
n=0のとき
$ [x^n]F = P(0)
これがループの終了条件
nが偶数の時
$ [x^n]F = [x^n]\mathcal{E}(F)
$ = [x^n](\frac{P}{Q} + \frac{P^*}{Q^*})/2
$ = [x^n](\frac{PQ^* + P^*Q}{QQ^*})/2
$ = [x^n](\frac{\mathcal{E}(PQ^*)}{QQ^*})/2
nが奇数の時
$ [x^n]F = [x^n]\mathcal{O}(F)
$ = [x^n](\frac{P}{Q} - \frac{P^*}{Q^*})/2
$ = [x^n](\frac{PQ^* - P^*Q}{QQ^*})/2
$ = [x^n](\frac{\mathcal{O}(PQ^*)}{QQ^*})/2
この時 $ QQ^*は偶関数
$ Q = q_0 + q_1x + q_2 x^2 + ...とする
$ Q^* = q_0 - q_1x + q_2 x^2 + ...を掛けると、奇数次の項が消えるから。
つまり数列としては、奇数番目が0の飛び飛びの列
これを圧縮したQ'が作れる
such as $ Q'(x^2) = QQ^*
形式的べき級数の積なのだから、数列の畳み込みで実装できて、飛び飛びで読めば良いだけ
同じように分子も飛び飛び数列なので圧縮できる
$ P_e(x^2) = \mathcal{E}(PQ^*)
$ P_o(x^2) = \mathcal{O}(PQ^*) / x
maspyさんはもっとコンパクトな分母を、おそらく手元での探索によって見つけて使っていた(勘違い)
が、オーパーツなので僕はここまでの理屈どおり$ \frac{x^5}{(1-x)^{16}(1+x)^5} を使った
下記のコードの[1,-2,0,2,-1]を見て勘違いしてたが、式は同じだった
code:python
den = np.array(1,-1, np.int64) for _ in range(5):
これは$ (1-x)^3(1+x) を展開したものだった
code:python
In 287: reduce(np.convolve, 1, -1 * 3 + 1, 1, 1) maspyさんはFFTを使った高速な畳み込みを使っていたけど、疲れたのでそちらを読むのはやめた
np.convoluteで済まそうとしたのだけど、int64をオーバーフローするので仕方なく自分で作った
Qm: $ Q^*
Qの奇数次の係数を符号反転するだけ
Qm = Q.copy() Qm[1::2] *= -1
code:python
def solve(N):
Q = reduce(
np.convolve,
1, -1 * 16 + 1, 1 * 5,
np.array(1, dtype=np.int64)) P = np.zeros(6, dtype=np.int64)
def conv(X, Y):
ret = np.zeros(n + m - 1, dtype=np.int64)
for i in range(n):
return ret
while N:
Qm = Q.copy()
QQm = conv(Q, Qm)
PQm = conv(P, Qm)
if N % 2:
else:
N //= 2